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(57) ABSTRACT 

Modeling and simulation of free and forced structural vibra- 
tions is essential to an overall structural health monitoring 
capability. In the various embodiments, a first principles 
finite-difference approach is adopted in modeling a structural 
subsystem such as a mechanical gear by solving elastody- 
namic equations in generalized curvilinear coordinates. Such 
a capability to generate a dynamic structural response is 
widely applicable in a variety of structural health monitoring 
systems. This capability (1) will lead to an understanding of 
the dynamic behavior of a structural system and hence its 
improved design, (2) will generate a sufficiently large space 
of normal and damage solutions that can be used by machine 
learning algorithms to detect anomalous system behavior and 
achieve a system design optimization and (3) will lead to an 
optimal sensor placement strategy, based on the identification 
of local stress maxima all over the domain. 

24 Claims, 28 Drawing Sheets 
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FINITE-DIFFERENCE SIMULATION AND 
VISUALIZATION OF ELASTODYNAMICS IN 
TIME-EVOLVING GENERALIZED 
CURVILINEAR COORDINATES 

RELATED APPLICATIONS 

This application claims priority to U.S. Provisional Patent 
Application No. 60/647,720, entitled “Finite-Difference 
Simulation and Visualization of Elastodynamics in Time- 
Evolving Generalized Curvilinear Coordinates,” by Upender 
Kaul, and filed on Jan. 1 9, 2005, which is hereby incorporated 
by reference in its entirety. The invention unifies an earlier 
invention, U.S. patent application Ser. No. 10/706,478, 
entitled “Enhanced Elliptic Grid Generation,” by Upender 
Kaul, and filed on Nov. 7, 2003, currently issued as U.S. Pat. 
No. 7,231,329 on May 23, 2007, which claims priority to U.S. 
Provisional Patent Application No. 60/425,750, entitled 
“Elimination of Parameter Input Requirement for Elliptic 
Grid Generation Methods in Engineering,” by Upender Kaul, 
and filed on Nov. 7, 2002 (both of which are hereby incorpo- 
rated by reference in their entirety), and a new solution 
method based on finite differences to simulate structural 
dynamic phenomena over time-varying grids in generalized 
curvilinear coordinates. 

STATEMENT OF GOVERNMENT INTEREST 

The invention described herein was made by an employee 
of the United States Government and may be manufactured 
and used by or for the Government of the United States of 
America for governmental purposes without payment of any 
royalties thereon or therefor. 

TECHNICAL FIELD OF THE INVENTION 

The present invention relates generally to numerical simu- 
lation of dynamic stresses, and in particular, to generalized 
curvilinear coordinate formulation for finite-difference pre- 
diction of stresses in elastic bodies under rotation. 

BACKGROUND OF THE INVENTION 

Conventional simulation methods to solve structural 
dynamics problems are in the domain of finite element tech- 
nology where the problem is solved in the modal domain and 
then the results are mapped into the time domain by appro- 
priate transformations. 

Limitations of the prior art have been in the difficulty in 
deriving new three-dimensional elements for different appli- 
cations of interest and the lack of ease in obtaining the tem- 
poral solution directly from the solution of governing elasto- 
dynamic pdes. 

The need to know the state of a structural system during its 
operation in terms of the physical output variables such as 
stresses and the geometric configuration of the system itself is 
essential for monitoring the system health. Such systems can 
be tested, prior to launching them in their operational domain, 
in a laboratory or through relatively inexpensive computa- 
tional simulations. 

Such systems when subjected to space and time varying 
loads during their operation can throw the system into unsafe 
states from the system’s health perspective. It is therefore 
essential to have a prior knowledge of such system states 
before the systems are commissioned. 

For the reasons stated above, and for other reasons stated 
below which will become apparent to those skilled in the art 
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upon reading and understanding the present specification, 
there is a need in the art for alternative modeling and simu- 
lation methods for structural systems. 

5 SUMMARY OF THE INVENTION 

Modeling and simulation of free and forced structural 
vibrations is an essential element of an overall health moni- 
toring capability for any structural system such as a rotorcraft 
io or any aerospace vehicle. In the various embodiments, a first 
principles finite-difference approach is adopted in modeling a 
structural subsystem such as a mechanical gear by solving 
elastodynamic equations in generalized curvilinear coordi- 
nates. While any structural subsystem can similarly be mod- 
1 5 eled, the various embodiments are described with reference to 
an annular disk, a thin solid disk and an idealized gear. Such 
a capability to generate a dynamic structural response has a 
wide applicability in a variety of structural health monitoring 
systems. Not only does this capability serve as a tool for 
20 understanding the dynamic behavior of a structural system 
and hence its improved design, but it also serves as a means by 
which a sufficiently large space of normal and damage solu- 
tions can be generated that can be used by a variety of 
machine learning algorithms to detect anomalous dynamic 
25 structural behavior of the system or to achieve a multi -func- 
tion design optimization of the given structural system. This 
capability will also aid in defining an optimal sensor place- 
ment configuration over structural subsystems for health 
monitoring, by identifying areas of local maxima of mechani- 
30 cal or thermal stress or loading. Such a capability to generate 
vibration response from a subsystem will also be useful in the 
area of vibration energy harvesting. Also, the methodology 
can be used to track stress wave propagation in a structural 
system which is useful in the health monitoring of such a 
35 system. 

The methodology is based on physics-based first prin- 
ciples, governing elastodynamics in the space-time domain. 
This innovation provides a powerful and yet simple method- 
ology to compute structural dynamic variables of interest 
40 such as stresses over an entire grid mapped over or inside a 
given body of interest directly in the time domain. The grid 
can be allowed to deform in time as the solution evolves. The 
simulation (deforming grids and stresses) can be visualized as 
the solution proceeds in time; the simulation can be sus- 
45 pended at any point in time based on the visualization of the 
state of the system and the simulation can be resumed or 
terminated altogether according as the evolving solution pro- 
ceeds within the expectation bounds dictated by physics. The 
attractiveness of the innovation lies in the intuitiveness of the 
50 approach where the physical variables such as stresses as well 
as the deforming body can be visualized directly in space and 
time, as the simulation proceeds. 

This simulation technology incorporates an innovative 
elliptic grid generation methodology that automatically 
55 updates the grid during the finite difference simulation of a 
given structural system directly in the time domain. The struc- 
tural simulation over such a geometry using elastodynamic 
partial differential equations (pdes) is itself innovative and 
gives results directly in the time domain. 

60 The various embodiments incorporate an enhanced elliptic 
grid generation algorithm (of the type described in U.S. 
patent application Ser. No. 10/706,478, currently issued as 
U.S. Pat. No. 7,231,329 on May 23, 2007) and a new three- 
dimensional finite-difference elastodynamic pde solution 
65 methodology in generalized coordinates. This software will 
enable simulation of dynamics of structural systems with 
deforming geometries directly in the time domain. An addi- 
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tional reference is “New Boundary Constraints for Elliptic 
Systems Used in Grid Generation Problems”, J. Computa- 
tional Physics, Vol. 189, 2003, pp 476-492, which is hereby 
incorporated by reference in its entirety. 

The solver of the various embodiments includes three main 
components. The first component is the enhanced elliptic grid 
generation algorithm that automatically updates the grid dur- 
ing simulation. The second component is the elastodynamic 
solver that solves the nine pdes, three for velocity components 
and six for stress components in generalized curvilinear coor- 
dinates over the grid generated by the first component. The 
third component is a visualization tool, such as the OPENGL 
based graphics and visualization software tool that animates 
the solution in time as the simulation proceeds. 

The various embodiments take any arbitrary three-dimen- 
sional geometry and pass it through the grid solver that gen- 
erates a smooth grid in generalized curvilinear coordinates 
about the geometry and then solve the nine elastodynamic 
pdes over each grid cell to yield the required structural solu- 
tion. Any updates, if needed, are made to the geometry, a new 
grid is generated about it, and the structural system is solved 
again. This process is continued for each time step over the 
required time of simulation. As the simulation proceeds, 
results are analyzed and if desired, the simulation is sus- 
pended and resumed or stopped altogether. Two novel and 
unique features of his invention include 1) the innovative 
enhanced elliptic grid generation algorithm under U.S. patent 
application Ser. No. 10/706,478; and 2) a new finite differ- 
ence based elastodynamic solver in three- -dimensional gener- 
alized curvilinear coordinates which solves for the structural 
solution directly in the time domain. The overall uniqueness 
of the innovation lies in that the methods of the various 
embodiments facilitate autonomous updating of the compu- 
tational grid in time as the simulation proceeds directly in the 
time domain. 

The various embodiments can be used to study a variety of 
structural and fluid- structural (in conjunction with a fluid flow 
solver) problems directly in the time domain in an intuitive 
fashion, with an ability to inspect and assess the simulation 
results as the simulation proceeds. There are a host of appli- 
cations for this innovation from space to aerospace to medical 
arenas, such as simulation of any structural system deployed 
in space or on earth and any artificial medical prosthesis in 
operation. 

For one embodiment, the invention provides a method of 
simulating time-dependent stress data for a structural system. 
The method includes generating an elliptic grid representa- 
tive of a geometry of the structural system, the grid compris- 
ing a plurality of grid cells, solving elastodynamic partial 
differential equations having velocity components and stress 
components in generalized curvilinear coordinates over each 
grid cell of the elliptic grid, updating the geometry of the 
structural system in response to solving the elastodynamic 
partial differential equations and generating a revised elliptic 
grid representative of the updated geometry of the structural 
system. For a further embodiment, the invention provides a 
computer-usable medium containing computer-readable 
instructions capable of causing a processor to perform the 
method. 

The invention still further provides methods and apparatus 
of varying scope. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. lisa representative grid of an annulus for use with an 
embodiment of the invention. 
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FIGS. 2A-2B are plots of radial and tangential stress for an 
annulus rotating at 100 rps comparing analytical and the 
computed results for two different grid resolutions in accor- 
dance with an embodiment of the invention at various grid 
5 dimensions. 

FIGS. 3 A-3B are plots of radial and tangential stress for an 
annulus rotating at 5 rps and 2000 rps, respectively, compar- 
ing analytical and the computed results in accordance with an 
embodiment of the invention at various speeds of rotation, 
to FIG. 4 is a representative grid of an annular disk for use 
with an embodiment of the invention. 

FIGS. 5A-5C are plots of radial and tangential stress for a 
thin annular disk rotating at 1 00 rps comparing analytical and 
the computed results at various axial stations in accordance 
15 with an embodiment of the invention. 

FIGS. 6A-6B are plots of stress for a thin solid disk com- 
paring analytical and the computed results in accordance with 
an embodiment of the invention at the quarter- thickness plane 
of a disk using various grid resolutions. 

20 FIGS. 7A-7B are plots of stress for a thin solid disk com- 
paring analytical and the computed results in accordance with 
an embodiment of the invention at the mid-plane of a disk 
using various grid resolutions. 

FIG. 8A is a cross-sectional view of a representative three- 
25 dimensional grid of a thin solid disk for use with an embodi- 
ment of the invention. 

FIG. 8B is a perspective view of the grid of FIG. 8A. 

FIG. 9 is a cross-sectional view of a representative three- 
dimensional grid of an idealized gear for use with an embodi- 
30 ment of the invention. 

FIGS. 10A-10B are plots of radial stress distributions 
along three “radial” lines of the idealized gear of FIG. 9. 

FIG. 11 is a cross-sectional view of a representative three- 
dimensional grid of an idealized gear for use with an embodi- 
35 ment of the invention and having a tooth region with reduced 
rigidity. 

FIGS. 12A-12C are plots of tangential stress signatures 
over time at various locations between the inner radius and 
outer radius at corresponding locations below a damaged 
40 tooth and a normal tooth of the idealized gear of FIG. 11. 

FIGS. 13A-13B are plots of radial stress signatures over 
time at various locations between the inner radius and outer 
radius at corresponding locations below a damaged tooth and 
a normal tooth of the idealized gear of FIG. 11. 

45 FIGS. 14A-14B are plots of shear stress signatures over 
time at various locations between the inner radius and outer 
radius at corresponding locations below a damaged tooth and 
a normal tooth of the idealized gear of FIG. 11. 

FIG. 15 is a plot of steady-state radial stress distribution of 
50 the idealized gear of FIG. 11. 

FIG. 16 is a plot showing a comparison between normal 
and damage shear stress vibration signatures midway 
between the inner and outer radii at tooth #10 of the idealized 
gearofFIG.il. 

55 

DETAILED DESCRIPTION OF THE INVENTION 

In the following detailed description of the present embodi- 
ments, reference is made to the accompanying drawings that 
60 form a part hereof, and in which is shown by way of illustra- 
tion specific embodiments in which the inventions may be 
practiced. These embodiments are described in sufficient 
detail to enable those skilled in the art to practice the inven- 
tion, and it is to be understood that other embodiments and 
65 variations may be utilized without departing from the scope 
of the present invention. The following detailed description 
is, therefore, not to be taken in a limiting sense, and the scope 
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of the present invention is defined only by the appended 
claims and equivalents thereof. 

Governing equation and boundary condition formulation 
in generalized curvilinear coordinates in three dimensions 
has been developed using the velocity- stress system of equa- 5 
tions of elastodynamics for bodies in rotation. The attendant 
partial differential equations are solved using a time-stag- 
gered leap-frog scheme. The stress boundary conditions in 
generalized coordinates are derived from a second order ten- 
sor transformation from the Cartesian coordinate system. The l o 
solution is based on first principles and does not involve 
lumped parameter or distributed parameter systems 
approach. The generalized curvilinear coordinate grids for 
the geometries considered here are created using a new and 
enhanced elliptic grid generation algorithm. The geometries 15 
considered herein as representative embodiments are an 
annulus, a thin annular disk, a thin solid disk and a thin 
multi-tooth gear. The solver, referred to herein as Finite- 
Difference Development of Linear Elasticity (FiDDLE), has 
been developed to solve the governing equations and generate 20 
the predictions. The predictions have been validated by com- 
paring them with the corresponding closed-form axisymmet- 
ric steady-state solutions for the annular and disk geometries. 
Then, predictions are made for a complex geometry of a 
multi-tooth thin spur gear in steady rotation as well as impul- 25 
sive rotation from an initial position of rest. The present 
methodology can be applied to study elastodynamics of com- 
plex shaped bodies under arbitrary dynamic loading. 

Nomenclature ^ 

E=Young’ s modulus of elasticity 
f=body force 

J=Jacobian of coordinate transformation 
k=bulk modulus 

M=metric coefficient matrix ^ 

R=right hand side vector 

q 7 , Q=velocity vector 

u,v,w=velocity components 

x=solution vector 

x,y,z, t=Cartesian coordinates and time coordinate 
^ 40 

5 77 =Kronecker delta 

X, (i=Lame constants 

V=gradient operator 

Q=rotational speed 

p=material density 45 

o^Poisson’s ratio 
T=symmetric stress tensor 

5, t|, £, e=generalized curvilinear coordinates and time 
coordinate 

Subscripts 50 

i, j, k=indices for coordinate directions 

Superscript 

T=vector transpose. 

The elliptic grid generation methodology makes it possible 55 
to generate two-dimensional and three-dimensional grids 
automatically around or inside arbitrary geometries, without 
any need for human intervention. For dynamically changing 
shapes, grids can be regenerated automatically during simu- 
lation as and when required. This would be pertinent to cases 60 
where a given geometry undergoes deformation such as bend- 
ing and twist, or where a crack propagates. In the present 
application of rotating gear component geometries, gear teeth 
could be subject to deformation, pitting, wear, and eventually 
cracks, and grids would need to be regenerated due to the 65 
corresponding changes in the boundary configuration. Hav- 
ing resolved the difficulty of automatically updating the grids 
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during a simulation, the finite-difference modeling of such 
problems has now become attractive. 

The need for numerical simulation of dynamic stresses 
over gears in mesh in both normal and damaged states has 
been delineated by the lack of any normal vibration data or 
any anomalous vibration data that may reflect the presence of 
gear damage such as in a pinion of an OH-58 helicopter 
transmission. Such simulated “clean” and “fault data” would 
aid in developing and enhancing fault-detection algorithms 
for such rotating systems by first, calibrating the damage- 
detection algorithms with the simulated clean vibration data 
and second, by validating these damage detection algorithms 
against simulated data corresponding to known damage. In 
these simulations, various forms of damage can be seeded and 
allowed to propagate in time and the corresponding data 
generated. In general, in such simulations, a vast variety of 
faults can be introduced into the system of interest, and the 
corresponding data could be recorded for use in the develop- 
mental work on fault-detection algorithms. 

Toward this end, as a first step, the present three-dimen- 
sional solver has been developed and validated by comparing 
its predictions with the known two-dimensional and three- 
dimensional theoretical steady-state solutions. The agree- 
ment has been shown to be good. 

Governing Equations 

The three-dimensional linear elastodynamic equations of 
motion describing the principle of momentum conservation 
and the constitutive equations governing the wave phenom- 
ena within an isotropic elastic body can be written as a system 
of nine equations, three for the velocities and six for the 
stresses, respectively. The velocity equations are given by: 

P dflrdjty+fi Eq. 1 

where the velocity vector, q z =(u, v, w). The body force 
vector is given by f.=(f x , f y , f z ), and the symmetric stress 
tensor, T y , has six distinct components. The stress tensor is 
expressed by the following tensorial equation: 

d t T iy =kd 0 - div Q+V-idfli+dfij) Eq. 2 

where div Q is the divergence of the velocity vector, Q (qj, 
h tJ is the Kronecker delta, and where X and p are the Lame 
constants — these elastic constants (e.g., \i=^z x Je xy is the rigid- 
ity or the elastic shear modulus, e xy is the strain caused by the 
stress x xy in the xy plane) characterizing the elastic behavior 
of the body are related to the Young’ s modulus of elasticity E, 
the Poisson’ s ratio a and the bulk modulus k by the following 
relations: 

is=ji(2fi+3X)/(fi+X) 

a=X/(2)(q+X)) 

k=E/(3(l-2o)) 

The elastic quantities ja and X are functions of space for a 
nonhomogeneous body. 

In the rotating frame of reference, the velocity equations 
(1) become: 

pdartyti+dii I Qxr\ 2 /2)+2ei Jk qjU> k +fi Eq. 3 

where IQI is a constant rotational speed (Q=c%), r=(x,y,z) 
is a positional vector, second and third terms on the right hand 
side represent the centrifugal and Coriolis forces associated 
with the rotating frame of reference. 

In generalized orthogonal curvilinear coordinates, Equa- 
tions 2 and 3, can be shown to assume the following flux- 
conservative form: 
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a.(A/J) 

+a 5 [(A/J)t,-(l/J)e, ( B+^C + t 2 D)] 

+a T ,[(A/J)n ( -(l/J)(ri J< B + ri J ,C+ri I D)] 

+ a,[(A/J)l r (l/J)^+% y C + % 7 D)]=R/J 
where 

A=[pU,pv,pw,T xx ,T^,T xz ,T >y ,T^,xJ r 
5 =/T xx ,T J(y ,T xz ,(X+2p)p,pv, i ww,X« J 0,X«/ r 
C=[x xy ,x yy ,x yz , Xv,/i/i,0 ,(X+2p)v^ w, Xv7 r 

D=ft xz ,r yz ,r zz , kw, 0 ,pp,^w Ji wv,(X+ 2 p) w/ r 

Jacobian of the transformation, J, is given by d(x,y,z)/d(|, 
t|£), and the metric quantities, ^ x , etc. have their usual 
meanings. The right hand side column vector, R, contains the 
centrifugal and Corioli s terms and additionally terms contain- 20 
ing spatial derivatives of the elastic constants, for nonhomo- 
geneous bodies. 

Physical quantities are normalized with the Young’s modu- 
lus, E, the acoustic speed, ^(X+2\i)/p, and the characteristic 
dimension such as the radius, r, for a solid disk or a shaft and 25 
(r OM ,-r /w ) for an annulus or an annular disk. 

Boundary Conditions 

The two-dimensional and three-dimensional validation 
examples presented here have corresponding axisymmetric 30 
steady-state theoretical solutions, and the predictions are 
directly compared with these theoretical solutions for a rotat- 
ing annulus (two-dimensional), rotating thin annular disk 
(three-dimensional) and a rotating thin solid disk (three-di- 
mensional). Boundary conditions corresponding to these 35 
closed-form solutions are transformed from the generalized 
coordinate space, (§,t|,Q, to Cartesian coordinates, (x,y,z), 
using contravariant tensor transformation. For example, sec- 
ond order stress tensor transformation between (x,y,z) space 
and (!,t),£) space is: 40 


5 

Eq. 4 
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dx l dx k 
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and the right-hand side vector, 

r=[ IVlI^IVlIlVrilT^JVlIlV^lT^IVril^, 

In a two-dimensional Cartesian coordinate system (y,z), 
the stresses would be related to their counterparts in the 
generalized coordinates (r|, £) as: 


= 2j[|V>;| 2 ^r m - 2|V,||V£tef 2 r„ f + \V Z\\\r cc ] 

r>z = 2 * hviji %t z r m + mwmyiz+v&yrx - ivtfw*] 

and 


J-zz = [|V Jj\ 2 ^yT m - 2|V,||V^^r^ + 1 Vf| 2 ^r ff ] 


where the Jacobian of the transformation, J', is given by 


d(x, y ) 


= 7!y£z-r}zty. 


and V is the gradient operator. Here, the orthogonal general- 
ized coordinates, (r|, £), correspond to the polar coordinates 
( 6 , r). 

In general, for a three-dimensional case, the stress tensors 
in the two coordinate systems would be related by the matrix 
form given above, i.e., x=M _ 1 R. 

The velocity boundary conditions are either of the 
Dirichlet type or of the Neumann type; the latter being 
derived from the governing equations for the velocity vector, 
once the stress tensor is updated at the boundaries according 
to the preceding formulation. 

Equation 4 in two-dimensional (y,z) or (r| , £) system can be 
written as: 

d^A =M'A^+L 'A^+R Eq. 4a 

where M’ and L 1 are matrices containing the variable p and 
various metric quantities, and R is the term containing body 
forces and any traction by way of boundary conditions. Equa- 
tion 4a can further be written as: 

dyA = UAU~ U^+L A-^+R Eq. 4b 


where x*=(£,'n,£), i, j span (§,r|,£) space and 1 , k span (x,y,z) 
space. The contravariant stress tensor transformation given 
above does not include normalizing factors. Using normal- 
ization, this transformation in matrix form can be written as: 

Mx=R 

where the metric coefficient matrix, 
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solution vector, 

X_ \yjocA’xyA'xz>X'xy?XyzA’zz\ 


or as: 

dyA=UX+L'A^+R Eq. 4c 

where X=AU - 1 A^ is a vector of two incoming character- 
istics, two outgoing characteristics and one neutral character- 
istic corresponding to the zero eigenvalue of the matrix M\ 
Therefore, on the boundaries, the solution variables are thus 

expressed in terms of these characteristics. 

55 

Numerical Method 

A second-order in time and space, time- staggered leap frog 
method is used to integrate the velocity and stress equations. 
Both the spatial derivatives and time derivatives are dis- 
60 cretized using central differences. A very small numerical 
damping term is used to eliminate the mesh drifting (check- 
erboard) instability. For one embodiment, a sampling rate of 
the order of 50 KHz is of interest, e.g., a 50 KHz sampling rate 
for the flight data of the OH-58 helicopter. For a further 
65 embodiment, the time step restriction imposed by the CFL 
condition of hyperbolic systems akin to explicit methods for 
integration of the dynamic system are of the same order as 
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those posed by the sampling rate of interest. For example, for 
a steel shaft of radius equal to 20 cm rotating at 100 revolu- 
tions per second (rps), a typical time step may be of the order 
of a microsecond, which translates to a simulation sampling 
rate on the order of 100 KHz. 

Results 

The predicted values of stresses shown are normalized by 
the Young’s modulus, E=2.1xl0 12 dyne/cm 2 , for the steel 
gear considered here. The lengths are normalized by the inner 
radius of the gear. Some steady-state theoretical results from 
linear elasticity were used as validation and verification base 
for the computations. Validation was done using three test 
cases: a rotating annulus (two-dimensional), a rotating thin 
annular disk (three-dimensional) and a rotating thin solid disk 
(three-dimensional). Then, predictions for a multi-tooth thin 
gear are presented. 

Rotating Annulus 

The thickness of the axisymmetric rotating annulus is 
assumed to be sufficiently small as compared to its radius so 
that the radial and tangential stresses do not vary over its 
thickness. This case is thus referred to as a rotating annulus 
here to avoid any confusion that may arise when we discuss 
the three-dimensional rotating annular disk with small finite 
thickness where the radial and tangential stresses do vary over 
the thickness of the disk. The closed form solution for the 
rotating annulus is given in Timoshenko, S. R and Goodier, J. 
N., Theory of Elasticity , McGraw-Hill, Inc., 1934. The 
weight of the annulus is neglected. FIG. 1 is a representative 
72x21 grid of an annulus for use with an embodiment of the 
invention. The annulus includes 72 grid points in the circum- 
ferential direction and 21 points in the radial direction. The 
annulus further has an inner and outer radii of 10 cm and 30 
cm, respectively. 

Results for an annulus of the type shown in FIG. 1 are 
shown for various grid dimensions and rotational speeds in 
FIGS. 2A-2B and 3A-3B. FIG. 2A shows a comparison 
between the analytical and the computed results at 100 rps. 
The abscissa shows the coordinate along any radial line from 
the inner radius to the outer radius, and the ordinate shows the 
normalized stress. The analytical, or theoretical, radial and 
tangential stress distributions along this radial ray are shown 
in solid and dotted lines respectively, and the computed radial 
and tangential stress distributions are shown in dashed and 
dotted lines with crosses respectively. Agreement between 
the computed and the analytical results is satisfactory. The 
grid used for this case is 72x11 , 1 1 points in the radial direc- 
tion and 72 points along the circumferential direction. As the 
grid is refined to 72x21, an improvement is seen in the pre- 
dictions, and the comparison with the theoretical results is 
good, as is shown in FIG. 2B. 

Results for the rotational speed of 5 rps with a grid of 
72x2 1 are shown in FIG. 3 A. Agreement between the predic- 
tions and the theoretical results is good, just as in the case of 
100 rps. Next, results for a high rotational speed of 2,000 rps 
with the same grid (72x21) are shown in FIG. 3B. Agreement 
between the predicted and theoretical radial stress distribu- 
tions continues to be good, but is subject to improvement with 
a finer grid. Grid dependence of the solution will be shown 
later in the case of a rotating thin solid disk. 

Rotating Thin Annular Disk 

An approximate theoretical solution of a rotating thin 
annular disk is given in Love, A. E. H., A Treatise on the 
Mathematical Theory of Elasticity, Dover, N.Y., 1944. The 
problem is treated as that of plane stress, in which the only 
nontrivial stress components are the radial and the tangential 
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stresses as in the case of the two-dimensional annulus, but 
now they also have a weak dependence on the thickness of the 
disk (along the axial direction). The thickness of the disk is 
taken to be 2.5 cm, while the inner and outer radii of the disk 
5 are 10 cm and 30 cm respectively, just as in the case of the 
two-dimensional annulus. The theoretical solutionis in defect 
near the ends, but is good at axial stations removed from the 
end planes. Here, the radial stress does not vanish along the 
outer or inner radii, as in the case of the annulus, but the 
to resultant radial tension between any two planes not too close 
to the end planes along the inner and outer radii vanishes. 
FIG. 4 is a representative 72x21x7 grid of an annular disk for 
use with an embodiment of the invention, with 7 points along 
the axial direction, 21 points along the radial direction and 72 
15 points along the circumferential direction. 

Results at three axial stations are shown in FIGS. 5A-5C. 
FIG. 5A shows the stresses at the plane next to the end plane; 
FIG. 5B shows the stresses at the section quarter thickness 
removed from the end plane; FIG. 5C shows the stresses at the 
20 mid- section of the disk. As before, predictions are shown with 
lines with crosses. The rotational speed of the disk is taken as 
100 rps, which is of the same order as the rpm rate of the 
pinion gear for an OH-58 helicopter transmission. The agree- 
ment between the predictions and the theory is quite good. 

Rotating Thin Solid Disk 

An approximate theoretical solution of a rotating solid thin 
disk is given in Timoshenko, S. R and Goodier, J. N., Theory 
of Elasticity, McGraw-Hill, Inc., 1934 and Love, A. E. H., A 
Treatise on the Mathematical Theory of Elasticity, Dover, 
N.Y., 1944. This case is similar to the thin annular disk with 
the difference that there is no inner radius here, and therefore 
the attendant vanishing radial tension condition at the inner 
radius becomes moot. 

35 A non-orthogonal cross-sectional grid is generated for the 
solid disk to avoid a polar coordinate singularity at the origin. 
Since the governing equations used in this study are strictly 
valid for an orthogonal curvilinear coordinate system, there is 
bound to be some discrepancy between the predictions and 
40 theory. The choice of the non-orthogonal cross-sectional grid 
in this case will thus help quantify the prediction errors asso- 
ciated with departure from a strictly orthogonal grid. 

The first grid, a coarse grid used for the thin solid disk is 
taken to be 21x21x7, with 7 points along the axial direction, 
45 21 points along the radial and circumferential directions . The 
radius of the disk is considered to be 10 cm and its thickness 
2.5 cm. The grids are progressively refined from 21 x21 x7 to 
41x41x11 to 61x61x11 to 81x81x11 to seethe improvement 
in the predictions. 

50 Boundary conditions at the end planes are imposed from 
theory. The rotational speed of the disk is taken as 100 rps. 

FIGS. 6A-6B show the predicted stresses at the quarter- 
thickness plane. FIG. 6A shows the radial stresses and FIG. 
6B shows the tangential stresses. Predictions are shown with 
55 lines with points and symbols and the theoretical stresses in 
solid line. Agreement between the predictions and the theo- 
retical results improves as the grid is refined from 2 1 x21 x7 to 
41x41x11 to 61x61x11 to 81x81x1 1. At 81x81x1 1 gridreso- 
lution, the predictions agree very well with the theory. As 
60 anticipated, any errors in the predictions at this stage may be 
due to the grid nonorthogonality, as discussed earlier, but it 
should be noted that the theoretical solution itself is also 
approximate. 

FIGS. 7A-7B show the predicted stresses at the mid-plane 
65 using the same grid resolutions as with FIGS. 6A-6B. FIG. 
7A shows the radial stresses and FIG. 7B shows the tangential 
stresses. Predictions are shown with lines with points and 
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symbols and the theoretical stresses in solid line. As expected, 
the predictions are progressively improved as the grid is 
refined. 

FIG. 8A is a 61x61 grid of a cross-sectional view of a 
representative 61x61x11 grid at any axial station of a thin 
solid disk for use with an embodiment of the invention. The 
three-dimensional grid for FIG. 8 A has a 61x61 grid and 1 1 
points along the axial direction. FIG. 8B is a perspective view 
of the grid of FIG. 8A. 

Rotating Thin Gear 

FIG. 9 is a cross-sectional grid of a 19-tooth gear for use 
with an embodiment of the invention. In the case of a gear, 
there is no theoretical solution available. Since the predic- 
tions have been validated for the annular geometries as dis- 
cussed above, results for a thin 19-tooth gear are presented 
here as a realistic application. There are some annulus solu- 
tion characteristics present in the gear predictions, as shown 
in the following figures. FIG. 10A shows a comparison 
among radial stress distributions along three “radial” lines, 
first one joining the inner radius with one end of the tooth at 
the base represented by a solid line, the second one between 
the inner radius and the top end of the tooth represented by a 
solid line with crosses, and the third one between the inner 
radius and the top middle point on the tooth. Similarly, FIG. 
10B shows the comparison among the three tangential stress 
distributions along these three “radial” lines. The stresses 
along the “radial” lines joining the inner radius and the top 
middle point on the tooth look similar to the results for the 
annulus. 

A physics based first principles approach is adopted to 
model and simulate vibration signatures from an idealized 
gear such as a thin spur gear of the type represented by FIG. 
9. The governing equations are solved using a finite-differ- 
ence approach as described above. The velocity- stress form 
of elastodynamic partial differential equations as used in 
earthquake signature modeling has been used here with the 
essential difference that in the present case, the velocity- 
stress system is solved in generalized curvilinear coordinates 
and the system being essentially dynamic entails the prescrip- 
tion of generalized characteristic boundary conditions based 
on the theory of hyperbolic systems; thus, artificial wave 
reflection and wave attenuation problems are totally elimi- 
nated. The vibration signatures are thus directly obtained in 
the time domain. A second-order accurate in time and space 
time-staggered leap-frog scheme is used to integrate the time- 
dependent partial differential equations. Idealized signatures, 
normal as well as damage vibration signatures, ensuing from 
an impulsive rotation of the gear are obtained and compared; 
normal signature is taken to be the one correspondent with 
homogeneous material properties throughout the domain. 
Damage signatures correspond to a case where the rigidity of 
one of the gear teeth is locally reduced. It is observed that 
significant deviations from the normal signature occur in 
amplitude and phase due to this damage. Using this approach, 
baseline or reference signatures can be obtained for any struc- 
tural subsystem which can be used to calibrate and validate 
various damage detection algorithms for such systems. 

Normal vibrations of an idealized multi -teeth steel gear are 
simulated by impulsively rotating the gear. This throws the 
gear into free vibrations about an equilibrium state that would 
be attained by it, if it had been set into rotation gradually from 
an initial state of rest. After about four rotations, the gear 
attains this steady state. 

The simulation considered corresponds to all the gear teeth 
except one, labeled tooth # 1 , having uniform material prop- 
erties as those of industrial steel. The shear modulus or the 
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rigidity of tooth #1 is decreased in a certain fashion over the 
region shown in FIG. 11 . This is just to mimic a damage state 
that would yield distinctly different vibration signatures from 
this particular tooth from those from the rest of the gear teeth. 
5 The gear is impulsively rotated at 6,000 rpm. The elastody- 
namic partial differential equations (pde), three for the veloc- 
ity vector and six for the symmetric stress tensor, are inte- 
grated in time, using fully characteristic boundary conditions 
as disclosed herein. An attractive element of the characteristic 
l o boundary condition approach is that the artificial wave attenu- 
ation and wave reflection problems associated with the tradi- 
tional boundary condition approach are entirely eliminated. 
The velocity-stress form of the elastodynamic pde has been 
used in geophysics to predict reference earthquake signa- 
ls tures. The integration is carried out until the end of the fourth 
rotation of the gear, when the vibrations have essentially died 
out and the equilibrium stress state is achieved. Thus the 
steady-state solution is obtained for radial, tangential and 
shear stress distribution all over the gear in constant rotation 
20 at 6,000 rpm. No grid independence study has been per- 
formed for this embodiment. But, having conducted the grid 
independence study using the present methodology for other 
geometries described above, the accuracy of the results pre- 
sented here is believed to be adequate in the present context. 
25 Also, as noted above, a small measure of Coriolis effect may 
be present even at 6,000 rpm during the dynamic state of the 
gear. 

The results from the simulation are compared side by side 
in terms of the time evolution of radial, tangential and shear 
30 stresses at selected locations on two selected gear teeth, tooth 
#1 and tooth #10, as shown in FIGS. 12A-12C, 13A-13B and 
14A-14B. Also shown is the steady state distribution of radial 
stress, all over the gear, with the damage tooth #1, in FIG. 15. 
The radial and tangential directions in the present case cor- 
35 respond to the two orthogonal generalized curvilinear coor- 
dinates, as shown in FIG. 11. 

A comparison of signatures from the damaged tooth #1 and 
a normal tooth #10 demonstrates that they vary significantly 
from each other. Also, certain insights can be drawn from 
40 simulations such as the present one that would aid in better 
design technologies of such systems. For example, as shown 
in FIG. 16, a comparison of time signatures of shear stresses 
corresponding to two cases, one where all the teeth are normal 
and the other where tooth #1 is damaged, reveals that the 
45 signal (shear stress wave) from the damaged tooth #1 reaches 
tooth #10 in about 0.3 millisecond, subsequent to an impul- 
sive rotation. 

These simulations can be used in conjunction with damage 
detection algorithms such as wavelets and machine learning 
50 methods for isolation, diagnosis and prognosis of system 
damage states. Such structural damage, amongst others, can 
be caused by unwanted material property variation due to 
thermal effects or manufacturing faults, sudden impact and 
repetitive loading, during the system operation. A wide vari- 
55 ety of these physics-based simulations could be carried out to 
compute various forms of normal and damage signatures. 
These signatures could then be used as reference signatures to 
calibrate and validate various damage detection algorithms 

60 CONCLUSION 

A new three-dimensional formulation in generalized cur- 
vilinear coordinates for the velocity-stress elastic system has 
been presented for a finite-difference solution over two and 
65 three-dimensional axisymmetric geometries as well as a gen- 
eral two-dimensional gear geometry. The axisymmetric pre- 
dictions are validated by comparing them with steady- state 
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axisymmetrical analytical solutions. Steady-state as well as 
unsteady predictions are also made for a thin multi-tooth gear. 
The various embodiments can be used to generate time-de- 
pendent stress data corresponding to a variety of boundary 
conditions on a variety of dynamical systems such as a rotat- 5 
ing gear simulating a gear pair in mesh, a rocket motor geom- 
etry and the space shuttle orbiter wing leading edge subjected 
to external dynamical loading. These numerically generated 
vibrations data can, in turn, be used to characterize different 
vibration states that the system may pass through. 10 

It is understood that the various embodiments may be 
performed by a processor in response to computer-readable 
instructions stored on a computer-usable medium, such as a 
fixed or removable storage media. Results may be presented 
to a user in any of a variety of formats, such as tabular, 15 
graphical or animation through the use of an appropriate user 
interface. Such computer systems are well understood. 

Although specific embodiments have been illustrated and 
described herein, it will be appreciated by those of ordinary 
skill in the art that any arrangement that is calculated to 20 
achieve the same purpose may be substituted for the specific 
embodiments described. Many adaptations of the invention 
will be apparent to those of ordinary skill in the art. Accord- 
ingly, this application is intended to cover any adaptations or 
variations of the invention. It is manifestly intended that this 25 
invention be limited only by the following claims and equiva- 
lents thereof. 

What is claimed is: 

1. A computer implemented method of simulating time- 
dependent stress data for a structural system, the method 30 
comprising: 

generating, using a computer system, an elliptic grid rep- 
resentative of a geometry of the structural system, the 
grid comprising a plurality of grid cells; 
solving, using the computer system, a three dimensional 35 
Finite Difference Model of elastodynamic partial differ- 
ential equations having velocity components and stress 
components in generalized curvilinear coordinates over 
each grid cell of the elliptic grid; 
updating, using the computer system, the geometry of the 40 
structural system in response to solving the elastody- 
namic partial differential equations; 
generating, using the computer system, a revised elliptic 
grid representative of the updated geometry of the struc- 
tural system; and 45 

animating, using the computer system, a solution of the 
elastodynamic partial differential equations in time as 
the simulation proceeds. 

2. The method of claim 1, wherein solving the elastody- 
namic partial differential equations further comprises solv- 50 
ing, using the computer system, the equations directly in the 
time domain. 

3. The method of claim 1, further comprising: 
introducing, using the computer system, a fault into the 

elliptic grid and allowing the fault to propagate in time 55 
while alternating solving the equations and generating 
revised elliptic grids. 

4. The method of claim 1, wherein the partial differential 

equations in generalized curvilinear coordinates, (!;, r| , €) in 

flux-conservative form are given by: 60 

dJA/T) 

dd(A/ J)1 - ( 1/ J) d^B+l^C+^D)] 

+<^ [(A/J)t| - ( l/J)(ri^B+r 1> ,C+r|^D)] 65 
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where R is a right hand side column vector containing cen- 
trifugal and Coriolis terms associated with the rotating frame 
of reference and, additionally, terms containing spatial 
derivatives of the elastic constants for non-homogeneous 
bodies, and J is the Jacobian of coordinate transformation, 
0(x,y,z)/3(^,r|,Q where d is the partial derivative with respect 
to the particular generalized coordinate variable r|, € € is 

the time coordinate in the generalized coordinate space and 
where the metric quantities ^ x , etc, have their usual mean- 
ings. The left hand side vectors are given by 

A=[pu, pv, pw, X xx , X yx , X xz , Xyy, X yz , Tj r 
B=[x xx , x xy , x xz , (X+2ji)u, pv, pw, Xu, 0, ku] T 
c=[x xy , Xyy, x yz , Xv, fju, 0, (X+2p)v, pw, Xv] r 

D=[x xz , x yz , x zz , Xw, Ofiu, Xw, }iv, (X+2p)w/ r 

where p is the material density of the structural medium; u, v, 
w are the three velocity components x tJ are the 6 components 
of the symmetric stress tensor; X and \i are the Lame con- 
stants. A is the solution vector, and B, C and D are the flux 
vectors 
Nomenclature 

J=Jacobian of coordinate transformation 
R=right hand side vector 
u,v,w=velocity components 

x,y,z, t=Cartesian coordinates and the time coordinate 
X, jii=Lame constants 
p=material density 
T=symmetric stress tensor 

T|, 5, e=generalized curvilinear coordinates and time 
coordinate 
Subscript 

i, j, k=indices for coordinate directions 
Superscript 
T=vector transpose. 

5. The method of claim 1, further comprising: 
generating, using the computer system, stress boundary 

conditions for the partial differential equations in gen- 
eralized curvilinear coordinates from a second order 
tensor transformation from a Cartesian coordinate sys- 
tem. 

6. The method of claim 1, wherein solving the elastody- 
namic partial differential equations further comprises inte- 
grating, using the computer system, the velocity and stress 
components of the equations using a second-order accurate in 
time and space, time-staggered leap frog method. 

7. The method of claim 1, further comprising: 
generating, using the computer system, vibration signa- 
tures from the time-dependent forcing functions. 

8. The method of claim 1, wherein the step of solving the 
three dimensional Finite Difference Model of elastodynamic 
partial differential equations includes solving, using the com- 
puter system, a set of nine elastodynamic partial differential 
equations, three partial differential equations for the three 
velocity components and six partial differential equations for 
the six stress tensor components. 

9. A computer implemented method of designing a struc- 
tural system comprising: 

generating, using a computer system, an elliptic grid rep- 
resentative of a geometry of the structural system, the 
grid comprising a plurality of grid cells; 
solving, using the computer system, a three dimensional 
Finite Difference Model of elastodynamic partial differ- 
ential equations having velocity components and stress 
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components in generalized curvilinear coordinates over 
each grid cell of the elliptic grid; 
identifying, using the computer system, areas of local 
maxima of stress directly at any given time in response to 
solving the elastodynamic partial differential equations; 
and 

modifying, using the computer system, the design of the 
structural system in response to identifying the areas of 
local maxima of stress. 

10. The method of claim 9, wherein modifying the design 
of the structural system in response to identifying the areas of 
local maxima of stress further comprises identifying, using 
the computer system, locations for sensor placement on the 
structural system corresponding to the areas of local maxima 
of stress for health monitoring of the structural system. 

11. The method of claim 9, further comprising: 
animating, using the computer system, a solution of the 

elastodynamic partial differential equations in time as 
the method proceeds. 

12. The method of claim 9, wherein solving the elastody- 
namic partial differential equations further comprises solv- 
ing, using the computer system, the equations directly in the 
time domain. 

13. The method of claim 9, wherein solving the elastody- 
namic partial differential equations further comprises inte- 
grating, using the computer system, the velocity and stress 
components of the equations using a second-order accurate in 
time and space, time-staggered leap frog method. 

14. The method of claim 9, wherein solving the three 
dimensional Finite Difference Model of elastodynamic par- 
tial differential equations includes solving, using the com- 
puter system, a set of nine elastodynamic partial differential 
equations, three elastodynamic partial differential equations 
for the three velocity components and six elastodynamic par- 
tial differential equations for the six stress tensor compo- 
nents. 

15. A computer implemented method of monitoring the 
health of a structural system, the method comprising: 

generating, using a computer system, an elliptic grid rep- 
resentative of a geometry of the structural system, the 
grid comprising a plurality of grid cells; 
solving, using the computer system, a three dimensional 
Finite Difference Model of elastodynamic partial differ- 
ential equations having velocity components and stress 
components in generalized curvilinear coordinates over 
each grid cell of the elliptic grid; 
generating, using the computer system, reference vibration 
signature corresponding to a given portion of the struc- 
tural system in response to solving the elastodynamic 
partial differential equations; and 
comparing, using the computer system, the reference 
vibration signature corresponding to the given portion of 
the structural system to a predicted vibration signature 
of the given portion of the structural system by a 
machine learning method; 

displaying, using the computer system, the reference vibra- 
tion signature and the measured vibration signature in 
real time during the operation of the structural system; 
and 

dynamically assessing, using the computer system, the 
health of the structural system by quantifying the devia- 
tion of the measured vibration signature from the refer- 
ence vibration signature directly in real time during the 
operation of the structural system. 


16. The method of claim 15, wherein solving the elastody- 
namic partial differential equations further comprises solv- 
ing, using the computer system, the equations directly in the 
time domain. 

5 17. The method of claim 15, further comprising: 

generating, using the computer system, stress boundary 
conditions for the partial differential equations in gen- 
eralized curvilinear coordinates from a second order 
tensor transformation from a Cartesian coordinate sys- 

10 tern. 

18. The method of claim 15, wherein solving the elastody- 
namic partial differential equations further comprises inte- 
grating, using the computer system, the velocity and stress 
components of the equations using a second-order accurate in 

15 time and space, time-staggered leap frog method. 

19. The method of claim 15, wherein solving the three 
dimensional Finite Difference Model of elastodynamic par- 
tial differential equations includes solving, using the com- 
puter system, a set of nine elastodynamic partial differential 

20 equations, three elastodynamic partial differential equations 
for the three velocity components and six elastodynamic par- 
tial differential equations for the six stress tensor compo- 
nents. 

20. The method of claim 15, wherein the reference vibra- 

25 tion signature is a reference displacement signature and the 

measured vibration signature is a measured displacement 
signature. 

21. The method of claim 15, wherein the reference vibra- 
tion signature is a reference acceleration signature and the 

30 measured vibration signature is a measured acceleration sig- 
nature. 

22. The method of claim 15, wherein the reference vibra- 
tion signature is a reference stress signature and the measured 
vibration signature is a measured stress signature. 

35 23 . A computer readable storage medium storing computer 

executable instructions which when executed on a computer 
to perform a method of simulating time-dependent stress data 
for a structural system, the medium comprising instructions 
for: 

40 generating an elliptic grid representative of a geometry of 
the structural system, the grid comprising a plurality of 
grid cells; 

solving a three dimensional Finite Difference Model of a 
set of nine elastodynamic partial differential equations 

45 having, three elastodynamic partial differential equa- 
tions for the three velocity components and six elasto- 
dynamic partial differential equations for the six stress 
tensor components in generalized curvilinear coordi- 
nates over each grid cell of the elliptic grid; 

50 updating the geometry of the structural system in response 

to solving the elastodynamic partial differential equa- 
tions; 

automatically generating in dynamic simulation time a 
revised elliptic grid representative of the updated geom- 

55 etry of the structural system; and 

animating, using the computer system, a solution of the 
elastodynamic partial differential equations in time as 
the simulation proceeds. 

24. A computer implemented method of designing a struc- 

60 tural health monitoring system that monitors the health of the 
structural system during its operation, comprising: 

generating, using a computer system, an elliptic grid rep- 
resentative of a geometry of the structural system, the 
grid comprising a plurality of grid cells; 

65 solving, using the computer system, a three dimensional 
Finite Difference Model of elastodynamic partial differ- 
ential equations having velocity components and stress 
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components in generalized curvilinear coordinates over 
each grid cell of the elliptic grid; 
identifying, using the computer system, areas of local 
maxima of stress directly at any given time in response to 
solving the elastodynamic partial differential equations ; 5 
identifying, using the computer system, locations for sen- 
sor placement corresponding to the areas of local 
maxima of stress for health monitoring of the structural 
system; and 

comparing, using the computer system, stress data from the to 
sensors with the three dimensional stress solution by 
using machine learning methodology; 


18 

displaying, using the computer system, the three dimen- 
sional stress solution and the stress data from the sensors 
in real time during the operation of the structural system; 
and 

determining, using the computer system, if operation of the 
structural system is safe or if the structural system is 
about to fail. 



